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Abstract. We use observations of the cosmological large-scale structure to derive 
limits on two-component hot dark matter consisting of mass-degenerate neutrinos and 
hadronic axions, both components having velocity dispersions corresponding to their 
respective decoupling temperatures. We restrict the data samples to the safely linear 
regime, in particular excluding the Lyman-a forest. Using standard Bayesian inference 
techniques we derive credible regions in the two-parameter space of nia and '^m^. 
Marginalising over '}2,my provides nia < 1.2 eV (95% C.L.). In the absence of axions 
the same data and methods give ^ < 0.65 eV (95% C.L.). We also derive limits on 
TOq for a range of axion-pion couplings up to one order of magnitude larger or smaller 
than the hadronic value. 
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1. Introduction 



The masses of the hghtest known particles (neutrinos) are best constrained by the largest 
known scales (the entire universe). The well-estabhshed method of using cosmological 
precision data to constrain the cosmic hot dark matter fraction [1, 2] has been extended 
to hypothetical low- mass particles, notably to axions, in several papers [3-5] . 

We return to this topic to extend previous studies by some of us [3, 4] in several 
ways. First, we update the cosmological data sets to include the Wilkinson Microwave 
Anisotropy Probe 3-year data as well as the baryon acoustic oscillations measurements 
from the Sloan Digital Sky Survey that have since become available. Second, we 
use standard Bayesian inference techniques to construct credible regions in parameter 
space, in contrast to the likelihood maximisation method used before [3,4]. Most 
importantly, we consider a two-component hot dark matter fraction consisting of axions 
and neutrinos. Since neutrinos are known to have nonvanishing masses, their hot- 
dark matter contribution is an unavoidable cosmological fit parameter. Axions and 
neutrinos decouple at different epochs and thus have different velocity dispersions that 
wc implement self-consistently. In this regard our work parallels a recent study by 
another group [5]. 

Wc begin in section 2 with a brief summary of the relevant axion parameters and 
their decoupling conditions. In section 3 wc describe the cosmological model and the 
parameter space we use. In section 4 we summarise the included data sets and briefly 
discuss our reasons for limiting the analysis to data in the safely linear regime of structure 
formation. We derive our new constraints in section 5 before concluding in section 6. 



2. Hot dcirk matter ctxions 



The Peccei-Quinn solution of the CP problem of strong interactions predicts the 
existence of axions, low-mass pseudoscalars that are very similar to neutral pions, except 
that their mass and interaction strengths are suppressed by a factor of order /jr/ /«, where 
/tt ~ 93 MeV is the pion decay constant, and /„ a large energy scale, the axion decay 
constant or Peccei-Quinn scale [6]. In more detail, the axion mass is 

^ _^ z'/' Um^ _^ 6.0 eV 

where z — mu/md is the mass ratio of the up and down quarks. We will follow the 
previous axion literature and assume a value z — 0.56 [7,8], but we note that it could 
vary in the range 0.3-0.6 [9]. Because of this uncertainty and to cover more general 
cases we will sometimes include a fudge factor Ca with the standard value 1. We will 
consider cases with —1 < log^Q^Ca) < +1. 

A large range of fa values (or, equivalently, rria values) can be excluded by 
experiments and by astrophysical and cosmological arguments [10]. Axions with 
a mass of order 10 /xeV could well be the cold dark matter of the universe [11] 
and if so will be found eventually by the ongoing ADMX experiment provided that 
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1 //eV < rria < 100 fieV [12]. In addition, a hot axion population is produced by thermal 
processes [13-15]. Axions attain thermal equilibrium at the QCD phase transition or 
later if /a < 10^ GeV, erasing the cold axion population produced earlier and providing 
a hot dark matter component instead. 

In principle, fa < 10^ GeV is excluded by the supernova SN 1987A neutrino burst 
duration [10]. However, the sparse data sample, our poor understanding of the nuclear 
medium in the supernova interior, and simple prudence suggest that one should not 
base far-reaching conclusions about the existence of axions in this parameter range on 
a single argument or experiment alone. Therefore, it remains important to tap other 
sources of information, especially if they are easily available. 

For those axion models with nonvanishing couplings to charged fermions, there exist 
stellar energy loss limits based on the axion-electron coupling that are competitive with 
the SN 1987A constraints so that here one does not rely on a single argument to exclude 
axions in the < 10^ GeV range. Therefore, we focus on hadronic models where axions 
do not directly couple to ordinary quarks and leptons. In this class of models all axion 
properties depend on /„ alone and not on model-dependent Peccei-Quinn charges of the 
ordinary quarks and leptons. 

If axions do not couple to charged leptons, the main thermalisation process in the 
post-QCD epoch is [13] 

a + n n + 7r . (2.2) 

The axion-pion interaction is given by a Lagrangian of the form [13] 

J^a. = (vrV+a^TT- + 7r°7r-a^7r+ - 27T+n-dy) d^a . (2.3) 

jTTja 

In hadronic axion models, the coupling constant is [13] 

Can = . ■ (2.4) 

"'^ 3(l + z) ^ ^ 

We note that in general the chiral symmetry-breaking Lagrangian gives rise to an 
additional piece for Can proportional to {rri^/ fT^fa) (7r°7r° -|- 27r~7r+) 7r°a. However, for 
hadronic axion models this term vanishes identically, in contrast, for example, to the 
DFSZ model (Roberto Peccei, private communication). 

Based on the axion-pion interaction, the axion decoupling temperature in the early 
universe was calculated by some of us in Ref. [4] , where all relevant details are reported. 
In our standard case we use the axion mass ma as our primary parameter from which 
we derive the corresponding axion-pion interaction strength by virtue of equations (2.1) 
and (2.2). Noting that even in hadronic axion models there is some uncertainty in this 
relationship due to the uncertain quark-mass ratio we consider also more general cases 
in which we include a fudge factor Ca as defined in equation (2.1), thus allowing for a 
more general relationship between rua and Can- 
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Table 1. Priors and standard values for the cosmological fit parameters considered in 
this work. All priors are uniform in the given intervals (i.e., top hat). 



Parameter 


Standard 


Prior 






0.01-0.99 






0.005-0.1 


h 




0.4-1.0 


T 




0.01-0.8 


ln(10i%) 




2.7 4.0 


ris 




0.5 1.5 


X^TOi. [eV] 





0-20 


rUa [eV] 





0-20 


logio(Ca) 





-1-1 



3. Cosmological model 

We consider a cosmological model with vanishing spatial curvature and adiabatic initial 
conditions, described by nine free parameters, 

^ = {a;dm,c^6,-f^o,T,ln(10^°A),ns,^m^,ma,logio(Ca)}. (3.1) 

Here, a;dm = ^dm^^ is the physical dark matter density, ujb — Q^/i^ the baryon density, 
Hq — h 100 km s~^ Mpc~^ the Hubble parameter, r the optical depth to reionisation, 
^4^ the amplitude of the primordial scalar power spectrum, and rig its spectral index. 
These six parameters represent the simplest parameter set necessary for a consistent 
interpretation of the currently available data. 

In addition, we allow for a nonzero sum of neutrino masses J^n^-u, a nonvanishing 
axion mass nia, and a fudge factor Ca relating nia to fa as defined in equation (2.1). 
These extra parameters will be varied one at a time, as well as in combination. Their 
"standard" values are given in table 1, along with the priors for all cosmological fit 
parameters considered here. 

4. Data 

4.1. Cosmic microwave background (CMB) 

We use CMB data from the Wilkinson Microwave Anisotropy Probe (WMAP) exper- 
iment after three years of observation [16-18]. The data analysis is performed using 
version 2 of the likelihood calculation package provided by the WMAP team on the 
LAMBDA homepage [19]. 

4-2. Large scale structure (LSS) 

We use the large-scale galaxy power spectra Pg{k) inferred from the luminous red galaxy 
(LRG) sample of the Sloan Digital Sky Survey (SDSS) [20, 21] and from the Two- 
degree Field Galaxy Redshift Survey (2dF) [22] . These power spectra are related to the 
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underlying matter power spectrum Pra{k) via Pg{k) = P{k)Pja{k), where the galaxy bias 
b{k) is conventionally assumed to be constant with respect to k over the scales probed 
by galaxy clustering surveys. 

However, recent studies suggest that this assumption may break down beyond 
k ^ 0.1 h Mpc"\ and may source the apparent tension between the SDSS and the 
2dF-inferred galaxy power spectra [20,23]. To model the effects of scale-dependent 
biasing when extracting cosmological parameters from galaxy clustering data, both the 
SDSS-LRG and the 2dF teams advocate the use of the Q^i fitting formula developed in 
Ref. [22] for ACDM cosmologies. See Ref. [24] for a detailed discussion. 

In the present work, however, we take the view that fitting formulae developed 
for standard cosmologies may not be applicable in nonstandard scenarios, particularly 
those involving new length scales arising from, e.g., axion and neutrino free-streaming. 
Developing an alternative formula to properly handle these nonstandard effects on 
the galaxy bias is also beyond our present scope. We therefore adopt a conservative 
approach, and use only power spectrum data well below A; ~ 0.1 /i Mpc^^, where a 
scale- independent bias is likely to hold true: 

• 2dF, A;inax ~ 0.09 h Mpc"^ (17 bands), 

• SDSS-LRG, k^^^ ~ 0.07 h Mpc'^ (11 bands). 

The combined set of these data is denoted LSS. We assume a scale-independent bias for 
each data set, and marginalise analytically over each bias parameter with a flat prior. 

4-3. Baryon acoustic oscillations (BAO) 

The baryon acoustic oscillations peak has been measured in the SDSS luminous red 
galaxy sample [25]. We use all 20 points in the two-point correlation data set 
supplied in Ref. [25] and the analysis procedure described therein, including power 
spectrum dewiggling, nonlinear corrections with the Halofit package [26], corrections 
for redshift-space distortion, and analytic marginalisation over the normalisation of the 
correlation function. 

4.4- Type la supernovae (SNIa) 

We use the luminosity distance measurements of distant type la supernovae provided 
by Davis et al. [27]. This sample is a compilation of supernovae measured by the 
Supernova Legacy Survey (SNLS) [28], the ESSENCE project [29], and the Hubble 
Space Telescope [30], as well as a set of 45 nearby supernovae. In total the sample 
contains 192 supernovae. 

4-5. Lyman-a forest (Lya) 

Measurements of the flux power spectrum of the Lyman-a forest has been used to 
reconstruct the matter power spectrum on small scales at large redshifts. By far the 
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Table 2. ID marginal 95% upper bounds on E and nia for several different choices 
of data sets and models. 



Data set 


Ca prior 


Em, [eV] 


rua [eV] 


WMAP+LSS+SNIa 


logio(a) = 


0.63 


2.0 


WMAP+LSS+SNIa+BAO 




0.59 


1.2 


Fixed E "ii^ = 








WMAP+LSS+SNIa+BAO 






1.4 


Fixed nia = 








WMAP+LSS+SNIa+BAO 




0.65 




WMAP+LSS+SNIa 


-1 < logio(a) < 1 


0.61 


2.2 


WMAP+LSS+SNIa+BAO 




0.60 


1.1 



largest sample of spectra comes from the SDSS survey. This data set was carefully 
analysed in McDonald et al. [31] and used to constrain the linear matter power 
spectrum. The derived linear fluctuation amplitude at A; = 0.009 km s~^ and z — 3 is 
= 0.452j:o:o6) and the effective spectral index ries = -2.32lj:o:o5- These results were 
derived using a very elaborate model of the local intergalactic medium in conjunction 
with hydrodynamic simulations. 

While the Lya data provide in principle a very powerful probe of the fluctuation 
amplitude on small scales, the question remains as to the level of systematic uncertainty 
in the result. The same data have been reanalysed by Seljak et al. [32] and Viel ct al. [33- 
35], with somewhat different results. Specifically, the normalisation found in Refs. [33- 
35] is lower than that reported in Ref. [31]. 

This question of normalisation is particularly important for bounds on the hot 
dark matter content of the universe. Since the free-streaming scale of light neutrinos 
or axions is larger than the length scale probed by Lya, their effect on the Lya data 
amounts to an overall change in the normalisation that is completely degenerate with 
any possible shift due to systematics. The Lja analysis in Ref. [31] already points 
to a higher fluctuation amplitude than that derived from the WMAP 3-year data; 
the addition of a hot dark matter component will render the two data sets even less 
compatible. This incompatibility in turn leads to a much stronger formal bound on the 
mass of the hot dark matter particle than would be expected considering the sensitivity 
of the present data (this is true for both neutrinos and other types of hot dark matter, 
such as axions). 

These considerations suggest that the Lya data are at present dominated by 
systematic effects. We therefore refrain from using them in the present analysis. 

5. Results 

We use standard Bayesian inference techniques, and explore the model parameter space 
with Monte Carlo Markov Chains (MCMC) generated using the pubhcly available 
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Figure 1. 2D marginal 68% and 95% contours in the ^ mi^-nia plane derived from 
the Ml data set WMAP+LSS+BAO+SNIa. The blue/solid lines correspond to the 
fudge factor being fixed at logio(Co) = 0, while the red/dashed hues indicate a top-hat 
prior on Ca in the interval —1 < log]^Q(Ca) < 1. 



CoSMoMC package [36, 37]. Our results are summarised in tables 2 and 3, and figures 1 
and 2. 

For the case of a standard hadronic axion our bounds on and are tabulated 
in table 2. When BAO data are included our bounds are almost identical to those 
recently derived in Ref. [5] based on their conservative data set. The main difference is 
that we do not use the HST prior on but instead include the SDSS-BAO data. Since 
the BAO data break the flm-h degeneracy, their inclusion has much the same effect as 
adding the HST prior. The importance of BAO data for the bound can be seen by the 
fact that the 95% upper bound is reduced from 2.0 eV to 1.2 eV. Our complete SNIa 
data set is also somewhat larger than the SNLS data set used in Ref. [5], containing in 
addition data from the GOODS and ESSENCE surveys. However, this has only a very 
modest impact on our results. 

Our neutrino mass bound < 0.65 eV (95% CL.) in the absence of axions 

is identical to that derived by some of us in Ref. [4], whereas the axion mass limit 
ma < 1.4 eV in the absence of neutrino masses found here is significantly weaker 
than the 1.05 eV limit found earlier [4]. The agreement of the neutrino mass limits is 
coincidental because here we use different data, notably excluding the Lyman-a forest, 
and a different statistical methodology (marginalisation instead of maximisation). The 
relative difference between the limits can be interpreted such that the axion bound 
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Figure 2. 2D marginal 68% and 95% contours in the 7na-\ogiQ{Ca) plane, assuming 
a top- hat prior on the fudge factor in the interval — 1 < log]^Q(Ca) < 1. 

benefits more from tlie inclusion of small-scale data than the neutrino mass bound, 
presumably because axions freeze out earlier and thus have a smaller velocity dispersion. 
Including the Lya data here would strongly improve both limits as can be gleaned from 
the results of Ref. [5]. Adding Lya to their conservative data set, the analysis of Ref. [5] 
finds that the marginalised axion mass limit improves by a factor 0.30, whereas the 
marginalised neutrino mass limit improves only by a factor 0.36, i.e., the relative gain 
for axions is 20% stronger. The changes in our new limits relative to those of Ref. [4] 
are in agreement with this picture. 

Returning to our new limits, an important observation is that the upper bound on 
the sum of neutrino masses is largely independent of whether or not massive axions are 
present. The 95% upper limit on X^^^v is in either case approximately 0.6 eV, a bound 
very close to that found in previous studies using roughly the same data combination 
[38-41]. 

In figure 2 we show how the bound on the axion mass changes as Ca is allowed 
to vary up or down by up to a factor of 10, assuming a uniform prior on \ogiQ{Ca) 
between —1 and +1. Note that the figure shows the 2D marginal contours, i.e., Ca and 
rua are fitted simultaneously. For \og^Q{Ca) < 0, the bound on does not depend on 
Ca because the number of degrees of freedom at decoupling, g^(To), is approximately 
constant for a large range of fa values (see table 2 of Ref. [4]). For logiQ^Ca) > 0, the 
value of g^{Tij) increases significantly with increasing Ca for a given m^. This increase 
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Table 3. ID marginal 68%/95% upper bounds on ^ rrii, and nia for fixed fudge factors 
Ca. The data set used is WMAP+LSS+SNIa+BAO. 







nia [cV] 


-1.0 


0.39/0.64 


0.51/0.98 


0.0 


0.37/0.59 


0.60/1.2 


1.0 


0.40/0.63 


0.69/1.4 



in g*{TD) leads to a corresponding drop in the present axion number density, 

Ua = ^^^^ X (5.1) 

For a fixed rUa this amounts to a decrease in the ratio Qa/Q^ with increasing Ca- The 
bound on rUa therefore becomes correspondingly weaker. 

As can be seen in figure 2, the 2D marginal 95% upper limit on rUa stays roughly 
constant at < 1.1 cV when log^Q(Ca) < 0, and increases roughly linearly with 
\ogiQ{Ca) to about 1.4 eV at \ogiQ{Ca) = 1. We stress again that the bounds in 
this figure are 2D bounds, and are formally — and often also in practice — not equivalent 
to ID bounds on rUa derived under the assumption of a fixed Ca- For instance, the 2D 
bound on at Ca = 1 in figure 2 is not exactly identical to the ID bound quoted 
in table 2 for a fixed Ca = 1. However, despite this formahty, the ma-logio(Ca) trend 
observed in figure 2 is also evident in table 3, which shows the ID marginal 68% and 
95% bounds on ma for fixed values of C^. 

Finally, we note that if Ca is increased much beyond 10, the bound will deteriorate 
rapidly because axion decoupling will have occurred beyond the QCD phase transition. 

6. Conclusions 

We have updated previous limits from cosmological structure formation on the mass of 
hot dark matter axions. This limit applies to axions which were thermalised, mainly by 
axion-pion interactions, in the early universe, and which subsequently decoupled from 
the thermal plasma while still relativistic. 

In the present study we investigate both the case where the neutrinos can be 
regarded as massless, as well as the case in which massive neutrinos are also allowed to 
contribute significantly to the hot dark matter fraction. In both cases we find an upper 
95% limit on the mass of hadronic axions of 1.1-1.2 eV when all available cosmological 
data, except the Lyman-a forest, are used. Reassuringly, we find that the bound on 
the sum of neutrino masses is almost completely unaffected by the presence of hot dark 
matter axions. 

Because of the uncertainty in the relation between the axion mass rUa and the energy 
scale /a, we have also studied the case in which the relation ma = 6.0 cV/ (/a/10^ GeV) 
is modified by a fudge factor Ca- We have studied Ca in the range 0.1-10, which is fairly 
representative of the model uncertainties. We find that the axion mass bound is largely 
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stable with respect to varying Ca in this range. For fixed values of C^, the ID marginal 
bound on goes from 0.98 eV at Ca = 0.1 to 1.4 cV at Ca = 10. Essentially, this 
means that for hadronic axion models the uncertainty of the light quark mass ratios 
have a negligible impact on the axion mass limit. 

Experimental and astrophysical limits on ma or /„ are always derived from limits 
on the axion coupling to different particles. The cosmological hot dark matter hmit, 
in contrast, primarily constrains the axion mass, with a very weak dependence on the 
axion-pion coupling. The hot dark matter limit of rua < 1 eV is very similar to the limit 
derived from globular cluster stars based on the axion-photon coupling. However, this 
coupling is quite uncertain even in hadronic models because even there it depends on the 
unknown electric charge of the heavy quark in KSVZ-type models. The hot dark matter 
limit implies that it is very difficult to escape the limit < 1 eV. One consequence is 
that in typical models, axions in the remaining allowed mass range necessarily escape 
freely from a supernova core. By courtesy of the SN 1987A neutrino burst duration, 
it follows that one can advance by another rung in the ladder of different limits and 
conclude that rUa < 10~^ eV [10]. While this SN 1987A energy- loss limit does not have 
an obvious loophole, we repeat that it is based on a very small sample of detected 
neutrinos and is subject to various nuclear-physics and axion-model uncertainties. 

Our results largely agree with those of Ref. [5] for their conservative data set. In 
contrast to Ref. [5] and to a previous study by some of us [4], wc have not included the 
Lyman-a forest data which could formally improve both the neutrino and axion mass 
limits roughly by a factor of 3. Wc have explained in section 4.5 that using the Lyman-a 
forest exposes one to the risk that large systematic uncertainties in the normalisation 
of the power spectrum at small scales may dominate the final result. 

The CAST experiment at CERN searches for axion-hke particles emitted by the 
Sun by virtue of their coupling to photons [42,43]. By including a helium fiUing of 
the magnet bores with variable pressure one can "adjust the photon mass," thereby 
allowing one to probe realistic combinations of and axion-photon coupling. The 
completed runs with '^He filling have already extended the experimental sensitivity to 
nia ~ 0.4 eV. Further extensions to up to rria ~ 1.16 eV with the forthcoming ^He runs 
over three years are on the agenda [44]. This search range is not excluded by our limits, 
particularly as we believe that more restrictive limits derived from the Lyman-a forest 
may be dominated by systematic effects that are not reliably controlled. 

One further caveat is that limits inferred from cosmological observations are by and 
large model-dependent. Additional free parameters not considered in this work, such as 
a nonstandard dark energy equation of state parameter, running in the primordial scalar 
spectral index, or a nonzero component of isocurvature modes in the initial conditions, 
could conceivably loosen the axion mass bound, as they have done many times before 
for the neutrino mass limit [38,41,45]. A significant and reliable improvement of 
cosmological hot dark matter limits is not immediately forthcoming. However, once 
data from the Planck CMB experiment [46] combined with other probes such as weak 
lensing surveys of galaxies [47-49] or of 21-cm emissions [50], or high-redshift galaxy 
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surveys [51, 52] become available, the sensitivity will be pushed down by as mush as an 
order of magnitude even in the face of more complicated cosmological model frameworks. 
In that event, a detection of axions by CAST in the vicinity of rUa ~ 1 eV will have 
important ramifications for observational cosmology. 
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